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ABSTRACT 


The structure of axi symmetric buoyancy-driven convection in a 

vertical cylinder heated from below is probed by finite element solution 
of the Boussinesq equations coupled with computed-implemented perturbation 
techniques for detecting and tracking multiple flows and for determining 
flow stability. Results are reported for fluids with Prandtl number of 
one and for cylinders with aspect ratio A (defined as the height to radius 
of the cylinder) between 0.5 and 2.25. Extensive calculations of the 
neutral stability curve for the static solution and of the nonlinear motions 
along the bifurcating flow families show a continuous evolution of the pri- 
mary cellular motion from a single toroidal cell to two and three cells 
nested radially in the cylinder, instead of the sharp transitions found for 
a cylinder with shear-free sidewalls. The smooth transitions in flow structure 
with Rayleigh number and A are explained by nonlinear connectivity between 
the first two bifurcating flow families formed either by a secondary bifurcation 
point for l\ < A* s 0.80 or by a limit point for A>_ A* • The transition 
between these two modes may be described by the theory of multiple limit point 
bifurcation. 
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1 . INTRODUCTION 


Since the early works of Benard (1901) and Rayleigh (1916) for a horizontal 
fluid layer, the onset and evolution of convective motions caused by temperature- 
induced buoyancy differences have been the focus of extensive theoretical and 
experimental research. For geometries with imposed temperature fields that are 
purely vertical and heated from below, convection begins at a critical tempera- 
ture difference, measured in terms of the Rayleigh number, beyond which the 
static fluid is unstable to small amplitude disturbances of the velocity, pressure 
and temperature fields. These critical Rayleigh numbers are determined as the 
eigenvalues for marginal stability in an analysis constructed from the Boussinesq 
equations linearized about the static state. 

The calculation of the fluid motions that evolve for Rayleigh numbers away 
from the critical values requires nonlinear analysis, either by perturbation 
methods (Schluter et al . 1965) or by numerical solution of the full Boussinesq 
equations. Perturbation methods are only feasible for geometries where at least 
part of the boundary is shear-free and eigenfunctions for the field variables 
can be written in terms of only a few special functions. Even in these systems 
limitations on the range of validity of the perturbation technique confine the 
results to flows only slightly perturbed from the rest state. Numerical calcula- 
tions have the potential for determining these flows over a much wider range of 
Rayleigh number and describing the nonlinear evolution of the structure with 
changes in the geometry of the cavity and Prandtl number. The purpose of this 
paper is to describe such a numerical study for axi symmetric convection in a 
vertical cylinder heated from below. 

We report calculations for the form and stability of the steady two- 
dimensional flows in a cylinder with rigid boundaries and insulated sidewall. 

These flows are computed by combining finite-e lenient methods for solving the 
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Bousslnesq equations with efficient computer-aided schemes for tracking the 
evolution and multiplicity of the flow field with changes in parameters, such 
as Rayleigh number. 

The numerical techniques for calculating solution multiplicity and 
stability have been developed as outgrowths of classical asymptotic analysis 
of bifurcation in systems of algebraic equations (cf. Keller 1977; Brown and 

Scriven, 1980a; Ungar and Brown 1982) and are applied here to the finite- 

/ 

dimensional equation set that results from the finite-element approximation 
to the Bousslnesq equations. With this approach, we are able to extend pre- 
vious calculations for the flows evolving from rest while simultaneously 
determining nonlinear interactions between families of flows. The stability 
results presented here are based on a linear analysis of the stability of the 
finite-element solutions to small perturbations in the field variables and 
are computed by schemes that make use of the connections between the change 
of stability of a family of flows and the occurence, of a critical value of 
Rayleigh where the flow is locally not unique. These techniques are general 
for all disturbances, except those leading to time-periodic bifurcation, and 
are much more efficient than the eigenvalue calculations used in previous 
studies (Brown and Scriven 1980a; 1980b). 

The cataloging of the branching and evolution of multiple flow fields is 
a fruitful approach to the description of nonlinear natural convection and the 
effect of varying parameters and geometry, as has been recently demonstrated 
by the studies of Daniels (1977 & 1979), Tavantis et al . (1978), and Hall and 
Walton (1977 & 1979). In this framework, the critical Rayleigh numbers {Ra^} 
determined from linear stability analysis mark the points of bifurcation between 
families of flows and the rest state; the eigenfunctions describe the structure 
of the new flows near the respective bifurcation points. A simple eigenvalue 
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Ra^ along the static family is then a simple bifurcation point joining 
two flow families, the static one and the new family of flows. These flows 
exist for some range of Rayleigh number beyond the Ra^ , as shown by asymp- 
totic analyses for rectangular cavities (Schluter et al . 1965; Daniel 1977; 

Hall and Walton 1977) and for a cylinder with shear-free boundaries (Liang 
et al. 1969; Rosenblat 1982). The critical values {Ra£*b have been obtained 
for vertical cylinders with several combinations of rigid and shear-free 
boundaries. 

Zierep (1959) first calculated these for a cylinder with rigid ends and 
shear-free sidewall and found that particular cellular flow structures were 
associated uniquely with each critical Rayleigh number. The ordering of the 
set {Ra^'b depended on the aspect ratio (defined as the radius over the 
height). Distinct values of A existed where two critical values of Ra were 
equal and where two flow patterns were equally likely. These values of Ray- 
leigh number are double bifurcation points (Iooss and Joseph 1980) in the full 
nonlinear problem and, as shown by Bauer et al. (1975), signal the possibility 
of secondary bifurcation along one of the new flow families for aspect ratios 
just slightly perturbed from this special value. Hall and Walton (1979) noted 
the existence of double points for the onset of convection in a rectangular 
cavity with rigid sidewalls and shear -free ends and used amplitude expansions 
to show the existence of a secondary bifurcation point. Double points are 
most simply found for a cylinder with shear-free boundaries all around (Liang 
et al_. 1969), as demonstrated by the plot of {Ra^b for axi symmetric modes of 
convection as a function of aspect ratio shown as Figure 1. Recently, Rosenblat 
(1982) has used approximations formed from truncated eigenfunction expansions to 
find secondary steady bifurcation of families of steady flows near double 
points between axi symmetric and non-axi symmetric modes for a cylinder with 
totally shear-free walls. 
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The boundary conditions on velocity at the sidewall of the cylinder 
have a striking effect on the evolution of the critical values with aspect 
ratio and thus on the existence of the double points and secondary bifurcation. 

As is shown below, a cylinder with rigid walls exhibits no crossing of the 
first and second critical Rayleigh numbers for aspect ratios between 0.5 and 
2.75. Instead, the cellular structure of the flow branching from Ra^ changes 
continuously with varying A. This result Is in agreement with the velocity 
profiles shown by Chari son and Sani (1970) for the onset of convection in the 
same problem and for calculations from an approximate analysis Brown and Stewart - 
son (1978) valid for large aspect ratios in a cylinder with rigid sidewall and 
shear-free ends. This latter linear analysis formed thebasis of nonlinear 
studies (Brown and Stewartson 1978 4 197S') of the effect on the structure of 
the bifurcating flow family of the sidewall and imperfections in the insulating 
condition along it. Exact calculations of the critical Rayleigh numbers for a 
cylinder with shear-free ends were reported by Joseph (1971), but insufficient 
information about the velocity component of the eigenfunction was given to deduce 
the evolution of the flow structure with A. 

Numerical calculations of the finite amplitude flows in a vertical cylinder 
have been performed with the three combinations of shear-free and no-slip boundary 
conditions mentioned above. Liang et, al . (1969) list sixteen finite difference 
calculations, all near the critical Rayleigh number, for different combinations 
of boundary conditions and including the dependence of viscosity on temperature. 
Jones et al . (1976) used time-dependent finite difference techniques to calculate 
the flow patterns in a cylinder with shear-free boundaries for Rayleigh numbers 
up to 100 Ra c and Prandtl numbers. ranging between 0.1 and the limiting case of 
infinity. They found a continuous family of flow fields that developed over this 
range of Ra with no qualitative change in the structure of the flow from the 
single toroidal cell predicted from the eigenfunction at Ra=Ra^. Chari son and 
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Sani (1975) have used approximations in terms of eigenfunctions to compute 
the finite amplitude flows and stability for a cylinder with all rigid 
boundaries, the configuration treated here. They present results for several 
aspect ratios and Prandtl numbers for Ra up to three times the critical value. 

The outline of the paper is as follows. The formulation of the natural 
convection problem Is presented in §2 and the finite element approximation 
and numerical techniques for tracking solution families are reviewed in §3. 

The links between bifurcations in the families of steady flows and the stability 
of each flow are established in §4 and the criteria for exchanges of stability 
are presented. Results of calculations of the onset of convection and the 
evolution of the flow families are represented in §5-7 . 
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2. FORMULATION 

We study the axisymmetric two-dimensional flows in a rigid vertical 
cylinder of height L and radius R filled with a fluid that has constant 
thermal diffusivity a and kinematic viscosity v. The two ends are taken 
to be isothermal with the lower one held at a temperature AT above the 
temperature T c of the upper surface. Two different thermal conditions are 
examined for the cylindrical sidewall; either it Is assumed to be a perfect 
insulator or it is taken as a perfect conductor so that a linear temperature 
profile exists along this wall connecting the temperatures at the top and 
bottom ends. Both conditions were studied in Chari son and Sani (1970 and 
1975} and the later case is included here for comparison to their calcula- 
tions. 

The dimensionless forms of the Boussinesq equations that govern the 
temperature e(r,z,t), pressure p(r,z,t) and velocity v_(r,z f t) fields are 

7 • v » 0 (2.1) 

f^- + v • Vv « -7p + Pr7 2 v + RaPr 6 e 2 , (2.2) 

-|| + v* ve * v 2 e , (2.3) 

where 7 is the gradient operator in cylindrical coordinates and e z is the 
unit vector in the vertical direction. The variables (y_,p,0) have been put 
in dimensionless form by scaling lengths with the height of the cylinder L, 
velocity with a/L, pressure with pa /L , time with L /a and temperature as 
(T(r,z ,t)-T c )/AT, where T(r,z,t) is the dimensional temperature field arid p 

t 

is the density of the fluid. The Rayleigh and Prandtl numbers appear in eqs. 
(2. 1-2.3) and are defined as 
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Ra ■ 0gL 3 AT/av , Pr s v/a , (2.4) 

where 3 and g are the coefficient of thermal expansion and acceleration of 
gravity, respectively.. 


The boundary conditions 

for the velocity field are 


v r ■ v 2 * 0 , 

0<r<_A , 

Z»0,1 

. (2.5) 

v r ■ v z * ° , 

r * A , 

0<z<l 

, (2.6) 

v r * 3v z /3r * 0 , 

o 

n 

L 

0<z<_l 

, (2.7) 


where AsR/L is the aspect ratio, When the sidewall is insulated the 
boundary conditions on temperature are 

6 * 0 0<r<A , z * 1 , (2.8) 

6 * 1 0<r<A , z * 0 , (2.9) 

c6/3r * 0 r * 0 , 0<z<l . (2.10) 

The thermal boundary conditions for a cylinder with a perfectly conducting 
sidewall are the same as eqs. (2.8-2.10), except the condition evaluated 
at r s A , eq. (2.10), is replaced by 

6 “ 1 -z , r * A , 0<z<l . (2.11) 

For either set of thermal conditions the entire equation set has the 
static solution 

y_ = 0 , 0 = l-z , p * p Q + RaPr(z-z 2 /2) , (2.12) 

from which the convective motions branch. We will represent these flows 
by Nusselt numbers defined as 



m 8 m 


Nu i 



( 2 . 13 ) 


where the temperature gradient Is evaluated along either the top (nJ) or 
bottom (N®) ends of the cylinder. 
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3, NUMERICAL METHODS 
3.1 Calculation of Steady Flows 

Finite element methods are well established for solving the steady and 
transient Bousslnesq equations and boundary condition*'. We use a method for 
mixed-order polynomial interpolation of velocities, temperature, and pressure 
that has been proposed before (see Zlenkiewlcz et al . 1976, Huyakorn et al. 

1979, Taylor and I jam 1979). The fluid domain (0<r<A, 0<z<1) is first divided 
into equally spaced quadrilateral elements. On this discretization, t ho com- 
ponents of velocity and the temperature are approximated by expansions in terms 
of biquadratic polynomials {<t> ^ (r ,z ) } as 

'v r (r,z, t)' 
v z (r,z,t) 

_e(r,z,t) 

(3.1) 

where the limits on the summations represent types of nodes in the elemental 
discretization; is the.number of nodes in the interior of the domain, 

N e represents the nodes along the top and bottom ends of the cylinder, and 

t , 

N s is the number of nodes along the sidewall. Each biquadratic function $ J (r,z) 
is defined to have a value of one at the node with which it is associated and 
to be zero at all other nodes. Because of this definition, the coefficients 
(u,(t), v . (t ) , 6 - ( t ) ) represent the values of the field variables at the loca- 

J J J 

tions of the nodes. The boundary conditions on temperature and velocity specify 
many of these coefficients and reduce the number of unknowns that must be compu- 
ted, with the exact number depending on the choice of thermal boundary condition. 
The remainder of the analysis is presented for the case of the sidewall being 
insulated. Then eqs. (2.8-2.10) set the coefficients {u®,v®, ®}, j*l , ...N e 
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and {u|,v|> , j»l , ... N s - 

The pressure is approximated by the expansion in terms of bilinear 
polynomials (t^(r,z)} as 

M 

p(r,z,t,) * l p i (t)y j (r,z) , (3.2) 

j*l 

where M is the number of vertex nodes in the finite element mesh. Each 
function ^(r,z) is defined so that it has a value of one at the jth 
vertex node and is zero at all others. More details of both sets of basis 
functions are available in many references (e.g, Thomasset 1981). 

The weak forms of the field equations are formed by applying Galerkin's 
method to eqs. (2.1 -2.3). Each field equation Is weighted with the basis 
function corresponding to the finite element expansion for the appropriate 
field variable and is integrated over the fluid domain: 


J i 

\ $ ve, dA * - 


' J 
V 


•Vv + £j<*Vp “ P^*v a v - RaPr6(e z *e k )3 dA 


, k*r,z , i=l , . . . N. 


(3.3) 


± 

at 


^ 6dA 

>v 


' * 

$ ] [v 2 6 - v*va ] dA 

h 


i*l .... N.+N 
i s 


(3.4) 




V V*v dA 


1 * 1 , 


M 


(3,5) 


The final forms of the residual equations are derived by applying the divergence 
theorem to eqs. (3.3, 3.4), by substituting in the expansions (3.1 , 3.2), and by 
incorporating the boundary conditions for temperature and velocity. After these 

t 

simplifications, the residual equations are reduced to the nonlinear ordinary 
differential equations 
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j*1 


1 du 
It 


au< 

I. -J- 


d'M dA * 


V 



9v 

r 

i/ 

3V 

r 

v r 

“Sr 

v z 

3z~ 

3& 1 * 

8v 

r 


3V 

r 

9r 

3 r 

3Z * 

~3z 1 


3v. 


~r. 2£. Pr izj:} 

9z Sr r Sr 1 


(3.6) 


N j dy ( 

l -3-ji. 4> * 4>^ dA * {$^(-v 

j*i dt b b 


9v_ 3V 3n , 3v„ 

’r Tr~ " v z Tz” “ sz “ Pr *F Tr " RaPr 


■ 3 * 3v z A 3* 1 3v * 


Pr( sFTr + 9r”# )}dA * 1-1, ... N i , (3.7) 


N i +N s de 


X dt 

j*l 


aA 


{^(-v — - v ZZ) 
' v r Sr v z sz ; 


S 6\ 




SO 1 99 iiL l&v ha 
9r Sr " Sz Sz* aM * 


M, ... N f +N s , (3.8) 


0 X 'i' j (Fi (r V ) 


!!= 

Sz 


-) dA = 0 , i=l 


(3.9) 


This set is compactly represented by the notation 


dx 


51 = R(x;Ra, Pr.A) 


(3.10) 


where x. is the vector of unknown coefficients 



' ” U N. ,V 1 ,v 2* 


* V Nj* 6 1 ’ e 2 ’ 6 N.+N c ,p l ,P 2’‘** P M^ » 

(3.11) 


and the vector R. represents the nonlinear algebraic equations that compose 
the right hand sides of eqs. (3. 5-3. 8). The mass matrix M is sparse and has 
components that are integrals of products of basis functions, but is mathema 
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ti cal 1 y singular because of the absence of time derivatives from the continuity 
equation. 

Steady state solutions of the natural convection problem are found by solv- 
ing the N t z (3N^ + N s + M) equations 

R(x;Ra) ■ 0 , (3.12) 


for given values of the parameters. We calculate families of solutions to eq. 
(3.11) by combining Newton's iteration with continuation methods (Kubieck 1976) 
for approximating the sensitivity of the solution vector x. on Rayleigh number. 
The details of these methods are laid out in Ungar and Brown (1982) and only 
points important to the developments in §4 are included here. 

Starting from a first approximation to the solution vector x_^ , Newton's 
method forms iterates as 


x< vfl > * x (i) -gj i+1 ^ s x (1 > - J" 1 (x* 1 ))R(x /’ J ; Ra) , (3.13) 


where J_ is the Jacobian matrix of the equation set (3.12), i.e. = aR^./BXj . 

The Jacobian matrix is asymmetric and has nonzeroes banded about the main diagonal 
because of the limited overlap of the finite element basis functions. The system 
of linear equations 

* -R(x (l) ;Ra) (3.14) 

is solved by Gaussian elimination with diagonal pivoting using the frontal routine 
of Hood (1976) to minimize the amount of computer memory necessary to perform the 
calculations. 

The first guess of the solution vector for a Rayleigh number R&ARa different 

t 

from Ra, where x_(R§) has been calculated before, is approximated as 
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x^'(Ra +ARa,Pr,A) * x(Ra Q ,Pr,,A) + Xfl a *<-Ra , 


(3.15) 


where the tangent vector x^ a is calculated from the linear set that results 
from taking the directional derivative of eqs. (3.12) along the solution 
family 


J(Ra Q ,Pr,A)x Ra * -(3R/3Ra) Ra 


(3.16) 


which is solved simultaneously with the last Newton iteration. 


3.2 Calculation of Limit Points, Bifurcation Points and Multiple Flow Fields 

Continuous families of flow fields are calculated by combining Newton 
iterations with continuation methods as long as the flows are uniquely speci- 
fied by Ra and no bifurcations to other families occur. The non-uniqueness 
with respect to Ra in a single family appears near a limit point Ra^ where 
the family reverses direction in Ra. Classical perturbation methods are well 
known for analyzing both limit and bifurcation points and have been adapted by 
several researchers (Keller 1977, Rheinboldt 1978, Brown et al . 1980) for im- 
plementation in numerical algorithms; the methods used here follow techniques 
developed by Keller (1977), 

Calculations near bifurcation and limit points are performed effectively 
by introducing an amplitude parameter e defined so that each family of flows 
is uniquely specified as (x(e), Ra(c)) in the neighborhood of the singular 
point fixed at e = 0 . Following perturbation methods used in classical bi- 
furcation analysis, the solution vector x and parameter Ra are expanded as 



9 


(3.17) 
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where the components (x^ , Ra^ are determined from the nonlinear equation 
set (3.12) expanded to the appropriate order of e . The sets governing the 
corrections (x^j , Ra (i)) and ^2)* Ra (2)^ are 


R x (x(0); Ra(0))2L (1 ) * -£* a (x(0h Ra(0))Ra (1) , (3.18) 

^(^(O); R&(0))x^ 2 ) * )£{i ) " 2 R x p a (x( 0 ); Ra(0) )'^ jRa^ j 

“B^ a (2L(°) ;Ra(0) )Ra^ 2 ) - EjiaRaRafi ) 

(3.19) 


In eq. (3.19), R xx is the (N t xN t xN t )-dimensional matrix with components 
{ R xx > i j k 5 3R i /Sx j 3x k ; R xRa is the (N t x N t )-d1mens1ona1 matrix with <R xRa >i j * 
SRVsx^sRa ; and R^^ is the vector with (B^ aRa > ^ 5 sV/SRa 2 . This last 
vector is zero for the Bousinesq equations. The (N^ x N^J-dimensional matrix 
R x appearing in eq. (3.18) and (3.19) is identical to the Jacobian matrix 
evaluated about the solution at the singular point. Terms in the expansion 
(3.17) can be computed ot.ce a suitable definition of e is selected; different 
definitions are used here for handling bifurcation and limit points. 

Riks (1972) first introduced and Keller (1977) refined the idea of intro- 
ducing the pseudo arc-length s for yielding local representations of the 
members of a solution family, even in the neighborhood of a limit point* and this 
is a suitable definition of the amplitude parameter e , that is e ■ s - s 0 . 

For computing around limit points, we introduce the arc-length evaluated from 
the known solution (x(s 0 ) , Ra(so)) through an additional residual equation 


r N 2 (s-s ) 2 +||x(s) t x_( s q ) 1 1 2 + |Ra(s)-Ra(s 0 )r * 


(3.20) 


’t+l 


2 u 2 

where ||x|^ is the R 2 -norm of x, i.e. ||x|g = J x^ . The tangent vector 
(x^, Ra^) for implementing continuation is computed by solving eq. (3.18) 
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augmented with the linearization of eq. (3.20) about the known solution 


(s-s 0 ) + d T x^i) + (Ra(s) - Ra(s 0 ) ) - 0 


(3.21 ) 


where s 3R N + -j/oX^ . The (N t +1 ) -dimensional matrix formed from eqs. 
t 

(3.18) and (3.21) is nonsingular at limit points (Keller 1977) and so the 
nonlinear set (3.12) augmented with eq. (3.20) is solved by Newton's method. 

The limiting value Rayleigh number Ra^ * Ra(s ]i ) is determined from the 
criterion Ra^j s (dRa/ds ) s s 0. At Ra * Ra^ the right hand side of eq. 

(3.18) vanishes and is determined as the solution of the homogeneous 
set and the orthogonal izati on condition (3.21); the solution of this problem 
is discussed in the next section. 

At values of Ra for bifurcation the Jacobian matrix is singular and the 
tanqent vector x^ a has multiple values. For simple bifurcation points only 
one eigenvalue of J_ passes through zero and bifurcation is detected by moni- 
toring the sign of the determinant of J^ . The null vector ^ corresponding to 
the zero eigenvalue and its adjoint vector y_ satisfy 

J(x(Ra c ); Ra c )z = 0 , i T (x(Ra c ); Ra c )y. = 0 , z T y.* 1 , (3.22) 

where J J is the transpose of J and is also banded. The null vectors are 
calculated by Gaussian elimination with partial pivoting as implemented in 
the frontal routine. The singularity of and the non-uniqueness of the null 
vectors are accounted for by ignoring the last row in the upper triangular 
matrix and by setting the N^-th element of z and £ to unity; these vectors 
are later normalized so that UzJ^ = jlyj^ s 1 • 

The continuation vector x^j written in terms of e , is decomposed at 
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Ra ■ Ra c into components in the null space and range of £ as 

2^1 j - Zq * Az_ , (3.23) 

where z is the null vector, A is an as-yet undetermined constant and z is the 

— ** U 

particular solution of eq, (3.18) that contains no component of z_ , i.e. that 
satisfies * 0 . Because the right hand side of (3.18) is homogeneous in 
, z^ s Ra^j£ where £ is the particular solution of 


i(x(Ra c ); Ra c )c = -R Ra (x(Ra c ); Ra c ) , (3.24) 

that is orthogonal to y. Equation (3.25) is solvable only if 

i. T SR a (x.(Ra ) c ; Ra c ) * 0 , (3.25) 


that is, the vector R^ a must be in the range of J. Bifurcation between the two 
solution families is guaranteed when Ra is a simple eigenvalue and eq. (3.25) 

V 

is satisfied (see Keller 1977 and Iooss and Joseph 1980). 

The tangent vectors to the base and bifurcating solution families are 
determined from (3.23) after A and Ra^ are computed from the definition of 
e and the solvability condition for the second-order problem (3.19). For analyz- 


ing simple bifurcation points in a flow family (x(e), u(e)) we use the amplitude 
definition 


s = Ay T (x_(e) - x(0)) + Ra^j(Ra(e) - Ra(0) ) , (3.26) 

where A and are associated with the solution family along which e is being 
defined. At first order in e eq. (3.26) requires that 



1 ■ A 2 + Ra^ j . 


(3.27) 
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The solvability condition for the second-order problem is 
A^C-j + ZARa^jCg + RapjC^ * 0 , 

where 

c i ■ *\x^ 

C 2 * l\ x £i + /S xRa I 

C 3 ' fc T W- + 

where the arguements of the derivatives of R. are the same as in (3.19). Equa- 
tions (3.27) and (3.28) are the same bifurcation equations derived by Keller 
(1977). The two pairs (A,.Ra^) are calculated as roots of these equations 
once the coefficients (C^> are determined at the bifurcation point Ra s Ra c . 

We approximate derivatives of the form R^ab. using the difference formula 

i\(x(0); Ra(0)ab = £ T [R x (x(0) + ia; Ra{0)) - R^xfO); Rs(0))]b/6 , (3.30) 

where o is a small positive constant, taken here as 10 . The matrix j^ Ra is 

calculated exactly. 

Flows in the bifurcating families are calculated by Newton iteration with 
an initial approximation given by eq. (3.17) . The amplitude s is picked to be 
large enough that the iterations converge to a flow in the new family. When 
Ra^j is zero the new family evolves either super- or sub-critically and we 
search values of Ra both above and below Ra c for new flows using first approxi- 
mations composed entirely of the null space. Either super- or sub-critical bi- 
furcation (Rap) = 0) occurs when C-pO •, see eq. (3.28), 
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4. BIFURCATION AND STABILITY 

The relative stability of flows in individual families to small perturba- 
tions of the velocity, temperature, and pressure fields can change only at bi- 
furcation and limit points in families of steady flows and at Hopf bifurcation 
points that mark the beginning of time-periodic motions. The discretized form 
of the Boussinesq equations (3.10) is a convenient framework for developing 
simple criteria for accessing the exchange of stability at the steady bifurca- 
tion and limit points detected by methods outlined in §3.2 . In this Section, 
we use ideas presented by Iooss and Joseph (1980) to extend earlier stability 
results developed by Szeto (1978). 

A 

The stability of a steady flow (x_(Ra), Ra) is analyzed with respect to 
small perturbations to the field variables represented in the same finite element 
bases presented in §3.1 as 

x_(t;Ra) = x(Ra) + u(t) = x(Ra) + ue at , (4.1) 

where the vector u, is independent of time. Linearizing eq. (3.10) for the 
evolution of a vector described by (4.1) gives the generalized eigenvalue problem 
for stability 

aMu_ * R v (x(Ra) ; Ra)u , (4.2) 

A 

where R x (x.(Ra); Ra) is the Jacobian matrix of R evaluated about the steady 

solution. The stability of this flow can be completely determined by computing 

the eigenvalues of eq. (4.2); however, repeated eigenvalue calculations are not 

feasible for the large matrices that result from the finite-element approximations. 

Instead, we develop stability criteria for the case when an eigenvalue a is zero 

at Ra * Ra and passes through zero as Ra is varied past Ra„ along the primary 
c c 

M 
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flow family. The criteria focus on determining the sign of the derivative 
(do/dRa) at Ra * Ra c . 

The adjoint eigenvalue problem to (4.2) is 

R^(x(Ra); Ra)y.* * cM T y* , (4.3) 

where is the adjoint eigenvector and a is the complex conjugate of a . 

At Ra * Ra , a is zero and the eigenvectors are the null vectors in eq. (3.22) 
with u. * z. and y* - y. . 

4.1 Exchange of Stability at a Limit Point 

Consider a family of flows parameterized by the amplitude e as (x(s), Ra(e)) 
near a limit point at Ra * Ra and corresponding to e * 0 . Expanding the 
eigenvalue and eigenvector in eq. (4.2) for small e, taking the inner-product 
with respect to y. , and evaluating at e * 0 yields 

fda-l . (A* Ra l}Ra n) , (4.4) 

e=0 y T Mz_ 

where the derivatives of R. are evaluated at x(0) and Ra(0). At a limit point, 
Ra^j is zero, x^ * Az. and (4.4) reduces to 

m . (4.5) 

LdeJ e*0 (i T Mz) 

The constant A is most efficiently eliminated from (4.5) by introducing eq. 

(3.26) as the definition of e which gives A s 1 . Equation (4.5) is simplified 
further by introducing the solvability condition for the second-order problem 
(3,19) evaluated at the limit point; 
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Substituting (4.6) into (4.5) gives the criterion for evaluating stability 
at a limit point as 



(4.7) 


We evaluate the sign ofRa^ from flows and their corresponding value 

' ' T T 

of e on both sides of the limit point and calculate the ratio I Br,/*. Hi 
using the adjoint null vector determined from eq. (3.22). The criterion 

(4.7) has been derived by Szeto (1978) and by Ungar and Brown (1983) using 
an alternative approach. 


4.2 Exchange of Stability at a Bifurcation Point 

Equation (4.4) gives the direction for the crossing of the eigenvalue 
through zero at either a bifurcation or limit point. Along the primary 
family, we assume that Ra ^ ^ is not zero and (4.4) is rewritten in terms 
of the coefficients in the bifurcation equation (3.28) 



AC.j + 82^2 / 1 \ 

(iV) 


(4.8) 


as first presented by Szeto (1978). The criterion (4.8) can be evaluated 
along either the primary or the bifurcating family using the pairs (A,Ra^) 
determined as roots of eqs. (3.27, 3.28) and is sufficient for determining 
exchanges of stability unless Ra^ is zero, that is, unless the bifurcation 
is super- or sub-critical. Then da/de * 0 along the bifurcating family and 
stability is determined by the second-derivative (d^a/de^) £= q . Szeto (1978) 
has developed a formula for this 'coefficient that involves the third derivative 
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-xxx eva ^ uate at the bifurcation point, a difficult quantity to calculate. 
Instead of using this result, we develop a more easily evaluated criterion 
based on the approach of Iooss and Joseph (1980) of linking the shapes of 
the bifurcating families to the exchange of stability along the primary 
family. 

2. Super- and Sub-Critical Bifurcation Points 

The equation set (3.10) is recast into local form expressed about the 
known, primary solution family (x(e), Ra(e)) by defining the new vector 

A 

w(e) * x(e) - x(e) and the equation set 

■ f(w; Ra(e)) a R(w + x(s); Ra(e)) - R(x(e); Ra(e)) , (4.9) 

Equations (4.9) have the steady solution w(e) * 0 and the vector f(w; Ra ( e) ) 
also satisfies the conditions 

-^4- (0; Ra(e) ) * 0 , (4.10) 

3 Ra* ” 

for all values of k. 

The form of the bifurcating family is analyzed by constructing an ex- 
pansion in e of w(e) and Ra(c) similar to eq. (3.17). In this local form, 
the first-order corrections are governed by 

f^wfe); Ra(e))Wg(e) + Ra e (e)f Ra (w(e) ; Ra(e)) s 0 , (4.11) 

with the subscript s standing for partial differentiation. When quantities 
in eq. (4.11) are evaluated at th'e bifurcation point (e = 0) this equation 
is analogous to eq. (3.18) . However, the condition eq. (4,10) with k s 1 
then reduces (4.11) to a homogeneous equation for w (0) 5 w^j(O) which has 
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a component only In the null space of the Jacobian matrix f (0; Ra(0)) s 

A A 

f x (x(0); Ra{0)) . Hence, the tangent vector in local form is, in some 
sense, strictly normal to the primary flow, that is 

w.(e) * z + 0(e) . (4.12) 

The stability of the flows along the bifurcating family (w(e), Ra(e)) is 
resolved in terms of an eigenvalue problem form by perturbing these steady 
solutions as w(t;e) * w(e) + £e' , where the magnitude of £ is small. 

This eigenvalue problem is 

f^wfeh Ra(e))£ * Y (s)Ml , (4.13) 

and its adjoint is 


fj"(w(e)i Ra(e))c* * y(e)M T £* » (4.14) 
where £* is ths adjoint eigenvector. 

An equation for the eigenvalue y(s) valid for any flow in the bifurcating 
family is derived by forming the inner product between and the first-order 
problem (4.11) and by using the relationship 


* y(e)|Vw € 

Then y(e) is given by 

, v ^‘loaWe); Ra(e)) 
y(e) = -Ra (e) T 

£ M ' > 


M £*w r 


(4.15) 


(4.16) 


A useful criterion for assessing stability at sub- or super-critical bifurcation 

t 

points is derived by expanding the right hand side of (4.16) for small e using 
the result (4.12) and the relationships Ra £ (e) s Ra^y + + 0(e 2 ) , 
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S? * £ + 0(e) » a nd f^ a <w(c) ; Ra(e)) » f wRa < 0 /, Ra(0)z e + 0(e 2 ) * 
R xRa (x(0); Ra(Q))ze+ 0(c 2 ) . Then eq, (4.16) reduces to 


Y(e) 


-Ra 


( 2 ) 


x \ Ra (i(°h s«(°))g 2 t 3) 

Jmz 


(4.17) 


which is re-written in terms of the slope of the eigenvalue a along the 
primary flow family as 


Y(e) 


■R 8 ( 2 ) (da/ dc ) e *o 2 

^7^T“ e 


+ 0(e 3 ) 


* 


(4.18) 


where (dRa/de ) e _q is the slope of the primary family at the bifurcation 
point and (da/de) eB0 is the slope of the critical eigenvalue. The result 
(4.18) is equivalent to the Factorization Theorem derived by Iooss and 
Joseph (1980) and gives the stability of the bifurcating family in terms 
of its evolution near the bifurcation point and the exchange of stability 
along the primary flow family. 
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5. CALCULATION OF THE ONSET OF CONVECTION 

The critical values of Rayleigh number where steady flows branch from 
the static state have been calculated by monitoring the determinant of the 
Jacobian matrix evaluated about this state, This determinant is plotted on 
Figure • 2 for a cylinder with A*0.5 and Insulated sidewall, as calculated with 
a mesh composed of four elements in each of the axial and radial directions. 

This mesh leads to a set of 224 nonlinear equations. The calculations were 
performed with Pr*1.0, however the critical values of Rayleigh number deter- 
mined are independent of Prandtl number. Although the determinant of ,J is 
extremely small in an absolute sense, its roots in terms of Ra are clearly 
seen in Figure 2 and can be systematically refined by numerical bisection. 

In most eases, we located the roots of det^) to within +5 in Rayleigh number. 

The lowest three critical values (Ra*^, Rap^, Rapp calculated here 
are compared In Table I to results reported by Chari son and Sani (1970) for 
cylinders with insulated sidewalls and aspect ratios between 0.5 and 2.70. 

The meshes used in our calculations, were varied so that the size of the elements 
remained essentially constant over this '"ange of aspect ratio. The finite 
element results agree well with those of Charlson and Sani, even for these 
relatively coarse discretizations. The structure of the flow field also is 
represented on Table I by the number of toroidal roll cells in the axial and 
radial directions; for example, 2R denotes a flow with two toroidal cells 
nested radially in the cylinder. For small A, our ordering of the second and 
third critical Rayleigh numbers does not agree with Charlson and Sani, who 
didn't report flows with more than one cell in the axial direction. 
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Increasing the number of elements in the mesh leads to more accurate 
results for the values of Ra,, and the’ corresponding eigenfunctions. This 
point is demonstrated by Table II where the lowest two critical values 
Ra^ are listed as a function of mesh for a cylinder with an in- 
sulated sidewall and A*1.0. Both values varied by less than 0.6% for meshes 
of between 4x4 and 8x8 elements. The values of Ra at the first several bi- 
furcation points are plotted on Figure 3 as a function of aspect ratio. 

The structure of the flow fields that evolve from these critical values 
are indicated by contours of the streamfunction shown for each curve and 
plotted more exactly on Figure 4. The flow field originating at the lowest 
value of Ra c evolves continuously with A from a single cell, to a two, and 
then to a three cell flow; see the plots for the lowest values of Ra c on 
Figure 4, 

For the calculations presented here with no-slip boundary conditions on 
all surfaces, there were no crossings of the curves of the lowest two values 
of Ra c ; the only sharp interchange of the ordering of two bifurcating families 
was found between the second and third critical values near AsA c =0.715 . The 
second critical point Ra£ 2 ) * Ra^) for this aspect ratio was a second-order 
critical point with two linearly independent eigenfunctions and hence had a 
Reisz Index of one (see Iooss and Joseph 1981). The role of this point in 
the nonlinear evolution of the flow structure is brought out in §7.2. 
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6. FINITE AMPLITUDE FLOW FIELDS: A-l .0 

The finite element approximations and computer-aided methods for 
tracking multiple solutions to algebraic equations have been used to 
calculate flow fields in the families emanating from the lowest four 
critical values of Rayleigh number for a cylinder with aspect ratio 
A=1 .0 and an insulated sidewall. Here and throughout the remainder 
of this report the fluid is taken to have a Prandtl number of one. Flow 
fields in each family close to Ra were found by using the first approx- 
imation eq. (3.23) with the corresponding eigenvector. Experience from 
weakly nonlinear analysis of Rayleigh-Benard problems (Schluter St al . ; 
Rosenblat 1982) and other finite amplitude calculations for the cylind- 
rical geometry (Liang et al . 1969; Charlson and Sani 1975) suggested that 
the bifurcating families evolve supercritical ly in Rayleigh number (towards 
increasing Ra). We searched for, and always found, the new flow families 
at (Ra-Ra ) greater than zero near Ra . The direction of the flows (either 

w w 

up or down along the centerline of the cylinder) in a particular family 
was arbitrary and pointed to two distinct families evolving from each bi- 
furcation point. Both were computed by changing the otherwise arbitrary 
sign of the constant A in eq. (3.23). 

The flow families computed with a 4x4 mesh are represented on Figure 
5 as a plot of the average Nusselt number between the top and bottom surfaces 
of the cylinder as a function of Ra. Since the sidewall was insulated, the 
heat fluxes through the top and bottom should have been the same and these 
two Nusselt numbers should have been equal. As discussed below, the discre- 

f 

tization error in the finite element approximation prevented this condition 
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from being satisfied. Flows composed of a single roll cell moving either 
up (HI) or down (ID) along the centerline developed for Ra>Ra c ; sample 
streamlines and isotherms from the ID family are shown in Figure 6. At 
large values of Ra-Ra. these flows developed a small secondary vortex in 
the upper corner of the cylinder which grew with increasing Ra. This 
secondary flow was found only for values of Rayleigh beyond those calcu- 
lated by Chari son and Sani (1975). 

The flow families that bifurcated from the second critical value 
Ra^ had two cells nested radially. Again there were two families of 
flows that differed only by the direction of the axial motion along the 

r 

centerline; we called these the 2RU and 2RD families, where the 2R desig- 
nation represented the radial structure of the flow. Sample streamfunctions 
and isotherms for flows in the 2RD family are displayed as Figure 7. As the 
Rayleigh number increased the purely radial orientation of the cells was 
lost and the flow evolved toward the same cellular configuration shown in 
Figure 6 for the ID family. 

Neither the (ID, 1U) or the (2RD, 2RU) families existed for Ra beyond 
a critical value Ra^*2.3xl0 4 , where the two families connected to form a con- 
tinuous solution curve; the value Ra=Ra^ is a limit point. The arguments 
for evaluating linear stability put forth in §4 show that the static solution 
lost stability at Ra a Ra^, as found by Charlson and Sani (1970), and that 
the 1U and ID flow families were stable. These flows lost stability 
at the limit point. The flows in the 2RU and 2RD families were all unstable. 
The stability of the flows in each f ami 1 i y is shown on Figure 5 by solid 
(stable) and dashed (unstable) curves. 

( 3 s ( a \ 

The flow families that bifurcated from the third Ra£ - and four Ra^ ' 
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critical values are also shown on Figure 5. Flows in the families evolving 

from Rap) had three radially nest cells (3RU and 3RD) and, passed quickly 

(41 

through a limit point in Rayleigh number. The flows associated with Ra^ ' 
had two roll cells stacked axially on top of each other. We called these 
families 2ASD and 2ASU, where the 2A designated the two axial cells, S stood 
for the plane of reflective symmetry through z=l/2 which divided the cells, 
and the D and U described whether the flow in the top cell was down or up 
along the centerline. Interestingly, these families were connected by a sec- 
ondary bifurcation point on the 2ASU family. Two other points of secondary 
bifurcation were also located on the 2ASD family, but because of the 2ASD 
flows were already unstable, the flow families evolving from these critical 
points were not calculated. 

The nonlinear connection of two families that appeared to be distinct in 
linear analysis was an important result of the results for A=1.0 and was found 
at all values of aspect ratio. We pause here to show that this phenomena was 
not an artifact of the coarse finite element approximation, but was indeed 
present for calculations with finer meshes. This point is made by examining 
Figure 8 where the results of tracking the ID and 2RD flow families are shown 
for four different finite element discretizations; clearly, all three sets of 
calculations were in qualitative agreement and the 6x6 and 8x8 meshes gave very 
similar results, indicating convergence of the calculation with mesh refinement. 

We also attempted to compare the Nusselt numbers calculated for flows in 
the ID family near the critcal value Ra^, 1 ) with those reported in Charlson and 
Sani (1975), but found the Nusselt numbers reported there to be much smaller 
than either the values we calculated or the values reported by Jones et al . 
(1976) for cylinders with large aspect ratio and shear-free sidewalls. Another 



comparison between our calculations and those of Chari son and Sani was made 
by calculating the flows that evolve from the lowest critical Rayleigh number 
for a cylinder with A*1.0 and a perfectly conducting sidewall . The Nusselt 
numbers evaluated at the top and bottom surfaces are plotted for this case 
in Figure 9 for two finite element grids along with the numbers reported by 
Charlson and Sani (1975). The sets of calculations are in good agreement, 
with the finite element results converging to those of Charlson and Sani for 
the finest mesh used. The ID family reached a limiting value of Ra for a 
cylinder with conducting sidewalls, just as it did for the perfectly insulated 
case; this case is discussed in more detail in another publication (Brown et-al. 
1983). 



7. EVOLUTION OF NONLINEAR FLOW STRUCTURE WITH ASPECT RATIO 


The double point found between the second and third critical Rayleigh num- 
bers at A=A C =0.715 (see Figure 3) hints that the nonlinear structure of the 
steady flows may undergo qualitative changes with varying aspect ratio. In 
this Section, we report results for aspect ratios between 0.5 and 2.00 for a 
fluid with Pr - 1.0 in cylinders with insulated sidewalls. The enormous number 
of calculations (almost 1500) needed to carry out this study has necessitated 
using the coarse finite element meshes listed in Table I for each value of 
A. Although calculations with these meshes may have errors in the overall 
heat balance of as much as ten percent for a few of the more vigorous flows, 
the checks of accuracy for A=1.0 discussed in §6 and mesh refinement of several 
of the calculations reported in this section give confidence that the quali- 
tative behavior of the flow fields are correct. 

7.1 The case A=0.5 

The flow families that evolved from the lowest two values of Ra„ are 

c 

represented on Figure 10 for A=0.5 and Pr = 1.0 . Just as for the cylinder 

with A=1 .0, flows composed of a single roll cell (1U and ID families) developed 

n ) t 

for Ra > Ra; ' up to a limiting value of Rayleigh numbers Ra s Ra . Repre- 

sentative streamlines for flows in the ID family are shown in Figure 11 and 

again developed a small secondary cell in the upper corner of the cylinder 

0 

which intensifies with increasing Ra up to Ra . The evolution of the ID and 
1U families past the limit point differed from the structure for the same flow 
families with A=1.0 discussed in §6. 

For A=0.5, the flow fields that bifurcated from Ra=Ra^^ had two cells 
stacked axially (2ASU and 2ASD) and also evolved toward higher values of Ra; 
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flows in the 2AUS family are shown in Figure 12. The plane of symmetry 
between the two cells was broken at secondary bifurcation points to new 
flow families, as indicated on Figure 11. The new flows branching from 
the 2ASU family are denoted as 2AUU and 2AUD, where the first U in both 
labels indicates that the top cell was circulating upward along the center- 
line and the suffices U and D indicate whether the top (U) or bottom (D) 
cell was the strongest. Members of the 2AUD family are displayed in Figure 
13 to show that velocity* fields away from the secondary bifurcation point 
had so large a bottom cell that the top one was pushed into the upper right 
hand corner of the cylinder and that the flows closely resembled those in 
the ID family. As implied by this remark, the ID and 2AUS families were con- 
nected along the solution curve that lead to the limit point Ra=Ra . 

The stability of the flows in these two families was determined entirely 

on the basis of the results for the stability of the flows emanating from 

0 

the static state, the numerical calculation of Ra as a simple limit point 
and the connectivity of the flow families. As shown on Figure 11, only the 
1U and ID flow families are stable in the sense described in §4. 

7.2 Flows near the double point: A=A C =0.715 

The double point at A* A c -0 .715 between the 2A and 2R families marked 
the first change in the flow structure from that shown for A* 0.5 . The flows 
computed for aspect ratios of 0.71 and 0.72, which bracket this double point, 
are represented in Figures 14 and 15. For A=0.71, the first two families have 
1U and 2A flow patterns and are connected at a secondary bifurcation point 
along the 2ASU branch, as was the case for A - 0.50 . The only change in the 
single-cell flows between A=0.5 and 0.71 was the development of a pair of 
limit points in this family and the creation of a new segment of stable flows; 



- 32 - 


this tehavior is depicted in Figure 14. At /fO .71 a secondary bifurcation 
point was found along the 2ASD family (two syirmetric axial cells with the 
top cell flowing down at the centerline) and the flows evolving from this 
point evolved continuously with increasing Rayleigh number into 2RU and 2RD 
f 1 ows . 

Changing the aspect ratio so that A>A C did not alter the connectivity 
of the 1U, 1R and 2ASU families, but switched the order of appearance of the 
2R and 2A families and moved the secondary bifurcation point from the 2ASD 
to the 2ASU branch. The change of the secondary bifurcation point as the 
apsect ratio passes through the double point A = A. is anticipated from the 

V 

asymptotic theories of Bauer et al . (1975) and Keener (1976) for bifurcation 
near such a second-order singularity. The connectivity between the 2R and 
2A flow families was not forecasted, although a similar structure occurred 
in the reaction-diffusion problem studied by Keener (1976). Figure 16 shows 
the evolution of the 2A and 2R families that we conjecture as A is varied 
through A . Because of the connectivity between the two families and the 
switching of the secondary bifurcation point from the 2ASD to the 2ASU fam- 
ilies only the 2AIJ and 2AD flows are thought to exist at A exactly equal to 
A c . 

7.3 Change in Flow Family Connectivity by Multiple Limit Point Bifurcation 

Varying the aspect ratio between 0.72 and 1.00 resulted in the loss of 
the secondary bifurcation point pictured in Figures 14-16 and the connectivity 
of the first two flow families shown in Figure 5. Evidence for the cause of 
this transition is displayed in the structure of the 2R and 2A flow families 

t 

for A=0.75 and 0.85 shown in Figure 17. At the lower aspect ratio the con- 
nectivity between the 2R and 2A was essentially the same as discussed for 
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A*0.72. However, at A*0 . 35 the 2R and 2A flow families were no longer con- 
nected and the two secondary bifurcation points along the 2ASU family were 
found to join a single family of flows. The 2R flows evolved continuously 
into the single-celled flows. 

The differences in structure between /f 0.75 and 0.85 suggests that the 

w 

flows emanating from the secondary bifurcation points interacted for A* A , 
0.75 <. A* £0.85 , in a way that replaced the coupling between the single- 
cell and 2A flows and resulted in the formation of continuous paths between 
the 1U, ID, 2RU and 2RD branches. Multiple-limit-point bifurcation, as dis- 
cussed by Decker and Keller (1980), is a possible mechanism for this transi- 
tion. The transitions through a multiple-limit-point are sketched in Figure 

★ 

18; here, the limit points in the 2R and 2AU families would coalesce at A* A 
and exchange connectivity, thereby leading to the structure shown for A>A C . 


7.4 The case A=2.0 

★ ... 

Increasing the aspect ratio beyond /?i\ did not change the connectivity 
of the flow families bifurcating from Ra ^ ^ and Ra^ shown f° r ^=0 ‘ ( see 

Figure 15). At large aspect ratio both families were composed of f s with 
two radially nested cells, which only differed in relative intense, stream- 
lines for sample flows in both families are shown in Figure 20. The flow 
patterns in two families that evolve from Ra£^ and Ra^ both have two radial 
cells but differ according to whether the inner or outermost cell is the most 
intense. Near the limit point, Ra=R e ~3.2xl0 4 , the two families converge with the 
outer cell driving the weaker inner cell to the centerline of the cylinder. 

A second set of flows bifurcating from Ra^ and Ra£ 4 ^ were alsu calcu- 
lated and evolved as a continuous loop with three radial cells (3R) in the 



- 34 - 


first family and four cells (4R) in the second. The coarse 10x4 mesh made 
it difficult to accurately resolve these velocity and temperature fields at 
Rayleigh numbers much above the critical values. The qualitative structure 
of the 3R and 4R flow families is most probably accurate; however, calcu- 
lations with more refined meshes are needed to accurately calculate the 
value for Ra at the limit point. 

Starting from an initial approximation and continuation vector calcu- 
lated for a flow belonging to the family that bifurcated from Ra^ , the 
Newton iterations converged to a flow in a new family. We tracked the 
family from this starting point and found it to evolve as a closed loop or 
isola which was not connected by bifurcation points to any other flows. Flows 
throughout the isola were composed of a dominant cell with motion up along 
the centerline, with a small counter-clockwise vortex in the lower corner of 
the cylinder. Sample streamlines along both the upper {with respect to 
Nusselt number) and lower branches of the isola are shown in Figure 21. The 
limit points which connected these two families mark changes in the relative 
stability of the two families, but the absolute stability of these flows could 
not be determined without either an understanding of the mechanism for creation 
of the isola from flow families originally connected to the static family or 
numerical calculation of the eigenvalues in eq. (4.2). 

We attempted to trace the evolution of the isola by numerically con- 

4 

tinuing solutions in the loop for A=2.0 at given Rayleigh numbers ( Ra=l . 0x1 0 
and 4xl0 4 ) to lower values of aspect ratio. An isola was successfully lo- 
cated for A s l .8 . Attempts to locate the isola at If 1.4 resulted only in 
calculation of flows in the 2R family. At this aspect ratio families of 
single-cell and 2R flows bifurcating from Ra^ and Ra,!^ were connected at 
a limit point RasRa^S.SxlO 4 , which this value for Ra^ was significantly 
higher than calculated for either A=1.0 or 2.0 . 
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4 4 

The flow fields calculated for Ra»l . 0x10 and 4.0x10 with 1.4<A<1.8 
are shown in Figure 22 and seemed to indicate that a continuous path in 
the solution space was traced. This result and the decrease in Ra^ between 
aspect ratios of 1.4 and 2.0 gave credence to the idea that the isola was 
created by the creation of singular points along the 2R family. The only 
candidate seemed to be the pinching of the loop to a point of sel ^intersection 
and finally to separation. This mechanism for isola generation has been docu- 
mented in models from reaction engineering (Uppal et al. 1976) and leads to 
the interesting conjecture that one branch of the isola may be composed of 
stable flows. We have not tried to analyze the creation of the isola beyond 
the results reported here. Again, finer finite element meshes may be necessary 
to resolve this detail. 
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8. DISCUSSION 

The most significant finding of this study was the connectivity pre- 
dicted between families of axisymmetric flows which originated at adjacent 
critical Rayleigh numbers. A continuous path of axisymmetric flows was 
found for all aspect ratios tested between 0.5 and 2.00 . The details 
of the path, i.e. whether it involved only a limit point or included sec- 
ondary bifurcation, did depend on aspect ratio. Imperfections in the ther- 
mal boundary conditions, which Introduce radial temperature gradients, rup- 
ture both the primary and secondary bifurcation points calculated here, as 
described by the analysis of Hall and Walton (1976), The connectivity be- 
tween the first and second axisymmetric families then is by a limit point 
along a continuous path for 0,5</,*2.70 . 

Connectivity between bifurcating families has not been detected in the 
previous asymptotic analyses of convection in either rectangular slits (Hall 
and Walton 1979) or cylinders with shear-free boundaries (Rosenblat 1982). 

Both analyses focused on values of aspect ratio for semi-simple double points. 
The first two double points between axisynmetric flows are represented on 
Figure 1 for a cylinder with shear-free boundaries. Neither of these second- 
order critical points exists for a cylinder with no-slip boundaries, as shown 
on Figure 3. We believe the imperfection in the spectrum of the linearized 
problem caused by varying the boundary conditions on the tangential velocity 
plays a major role in establishing the connectivity between solution families 
observed here, but not seen in calculations for cylinders with shear-free 
sidewalls (Jones et al. 1976; Brown 1983). A systematic asymptotic study of 
convection in cylinders with slightly sticky boundaries is underway. 

The calculations presented where are only the first step toward a compre- 
hensive understanding of convective transitions in a cylinder heated from below. 
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Three-dimensional and time-periodic flows are almost unexplored theoretically. 
Nonaxi symmetric convection modes are known (Chari son and Sani 1971) to be most 
dangerous at large and small values of A of a rigid cylinder. Even when an 
axlsymmetric mode is the first to bifurcate, as is the case for Pr*1.0 and 
A*1 . 0 , three-dimensional convection may dominate after a secondary bifurcation, 
as demonstrated by a few calculations in Charlson and Sani (1975) for a cylinder 
with conducting sidewalls. Three-dimensional convection patterns and time per- 
iodic flows have been reported by Olson and Rosenberger (1979) for gases in a 
cylinder with A*l/6 heated from below for Rayleigh numbers 5.86 times Rajp^ . 

The numerical value of Ra for the onset of oscillations is no doubt a function 
of both aspect ratio and Prandtl number. 

The strategy presented here for computing steady flows and analyzing non- 
linear structure and stability generalizes to the study of three-dimensional 
flows by the finite element method and to other numerical approximations for 
computing two-dimensional convection which incorporate Newton's method for 
solution of the full set of residual equations (van Steeg and Wesseling 1978; 
McDonough and Catton 1982). Only the availability of large, fast super - 
computers limits this approach, as it did for the pioneering calculations of 
Charlson and Sani (1975). Nonaxi symmetric calculations are now feasible in 
geometries where the boundary conditions allow spectral representation of the 
azimuthal dependence of the field variables. New numerical approximations that 
are a hybrid of spectral and finite element approximations are being developed 
for calculating steady three-dimensional convection in a cylinder. 

Just as in asymptotic expansions based on eigen-modes, care must be taken 
to access the range of validity, of the numerical appoximation at large values 
of jRa-Ra | or when the flow pattern involves multiple cells. Spurious steady 
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solutions to discretized equations have been reported in several studies 
(Schri eber and Keller 1983} Chang and Brown 1983} Yeh et al. 1983) when 
the numerical approximation, either finite difference or finite element, 
cannot resolve boundary layers and separation in a flow field. Similar 
numerical artifacts can result at low values of |Ra-Ra c | when the finite 
element approximation is insufficient to resolve multiple flow cells} this 
limitation has restricted our calculations to A less than 3. Spectral de- 
compositions are better suited to calculations at large A. 

The authors are grateful to the Materials Processing Program of the 
U.S. National Aeronautics and Space Administration and to the Information 
Processing Services at Massachusetts Institute of Technology for support 
of this work. Y. Yamaguchi was supported by Mitsubishi Chemical Industries 
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FIGURE CAPTIONS 


Table I. 


Figure 1. 


Figure 2. 


Figure 3. 


Figure 4. 


Figure 5. 


Figure 6. 


Figure 7. 


Comparison of first three critical Rayleigh numbers calculated 
by finite element analysis with values reported by Charlson and 
Sani (1970) for cylinders with aspect ratios between 0.5 and 
2.75 and insulated sidewalls. The mesh was adjusted to conform 
with the aspect ratio of the cylinder. 

Schematic of critical Rayleigh numbers (Ra^> as a function of 
aspect ratio A for a cylinder with shear-free boundaries. Each 
curve is associated with a particular flow structure and the 
ordering of the modes is interchanged by crossing the critical 
values A q 

Determinant of the Jacobian matrix evaluated about the static 
state as a function of Rayleigh number for A=0.5 and for a 4x4 
finite element mesh. The cylinder has insulated sidewalls. 

Lowest several critical Rayleigh numbers for a cylinder with in- 
sulated sidewalls as a function of aspect ratio A. 

Contours of streamf unction for the eigenfunction corresponding 
to the lowest three critical values of {RaJ^H for aspect ratios 

V 

between 0,5 and 2.75. 

Families of axisymmetric flow fields in a cylinder with insulated 
sidewalls for A=1.0 and Pr=1.0 . Stable flows are shown by solid 
( — ) curves and unstable flows by dashed (---) ones. 
Representative streamlines and isotherms for flows in the ID 
family that occur before, the limiting value of Ra for AfI.O . 
Representative streamlines and isotherms for flows in the 2RD 
family for A«1.0 . Each of these flows are unstable. 


Figure 8. 


Figure 9. 


Figure 10. 


Figure 11 . 
Figure 12. 
Figure 13. 


Figure 14. 


Figure 15. 
Figure 16. 
Figure 17. 


Families of ID and 2RD flows in a cylinder with insulated side- 
walls computed with four different finite element girds; A*1.0 
and Pr*l .0 . 

Nusselt numbers at top and bottom surfaces of a cylinder with 
perfectly conducting sidewall; A*1.0 and Pr*l»0 . Finite ele- 
ment results for two meshes are shown along with the results of 
Charlson and Sant (1975). 

Families of axisymmetric flow fields in a cylinder with insulated 
sidewalls for A a 0.5 and Pr=1.0 . Stable flows are denoted by 
solid ( — ) curves and unstable flows by dashed (---) ones. 
Representative streamlines and isotherms for flows in the ID 
family that occur before the limiting value Ra^ . 

Representative streamlines and isotherms for flows in the 2ASU 
family. 

Streamlines and isotherms for members of the 2AUD family which 
show the evolution of the flow field with changing Ra into the 
ID family; compare with Figure 11. 

Families of axisymmetric flow fields in a cylinder with insulated 
sidewalls for A=0.71 and Pr--1.0 . Stable flows are denoted by 

solid ( ) curves and unstable flows by dashed ( — ) ones. 

Families of axisymmetric flow fields in a cylinder with insulated 
sidewalls for A=0.72 and Pr=1.0 . 

Evolution of flow families between families with two axial and 
two radial cells as aspect ratio is varied through A c - 
Structure of 2ASD, 2ASU'and 2R flow families for A-0.75 and 0.85 


Pr a 1 .0 



Figure 18. 


Figure 19. 
Figure 20. 


Figure 21 . 


Figure 22. 


Evolution of multi pi e-limit point bifurcation with varying aspect 
ratio through A s Af which is proposed as the mechanism for separation 
of the 2R and 2A flow families. 

Families of axi symmetric flow fields in a cylinder with insulated 
sidewalls for A =2.0 and Pr*1.0 . 

Representative streamlines for flows in the 2RU families for A*2.Q 
and Pr*1.0 ; the six sets of contours (a-f) follow the evolution of 
the flow from Ra£^ around the limit point at RasRa^-S^xlO^ . 
Representative streamlines for flows in the isola found for A *2.0 
and Pr=1.0 ; the six sets of contours (a-f) follow the evolution of 
the flow from the branch with higher values of average Nusselt num- 
ber to the lower branch. 

Streamlines of flows tracked for Pr*l .0 and Rayleigh numbers of 
4 4 

1.0x10 and 4.0x10 with changing aspect ratio. 
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